In [1]:
import re
import glob
from IPython.core.display import Image
In [2]:
Image(filename='figures/experiment-setup.png', width="100%")
Out[2]:
Rohde & Schwarz SMBV100A vector signal generator,
Mini-Circuits 50 ohm 30 dB attenuator,
coaxial cable
In the following, each method is marked with an abbreviation as follows:
Three settings for the USRP were chosen:
$t_s$ = 25.0 ms, $f_s$ = 1 MHz, $N_s$ = 25 kSamples
Low-end embedded profile. Limited by the ADC performance and available RAM for ADC buffer. This is expected to be the maximum performance using single 1 MHz ADC on a STM32 microcontroller with 100 kB of RAM.
$t_s$ = 12.5 ms, $f_s$ = 2 MHz, $N_s$ = 25 kSamples
High-end embedded profile. Limited by the ADC performance and available RAM for ADC buffer. This is expected to be the maximum performance using dual 1 MHz ADC in interleaved mode on a STM32 microcontroller with 100 kB of RAM.
$t_s$ = 10.0 ms, $f_s$ = 10 MHz, $N_s$ = 100 kSamples
Profile for a high-performance CPU and SDR front-end.
SNE-ISMTV contains an analog energy detector. Sampling rate is fixed at $f_{sSNE}$.
In [3]:
fsSNE = 1/6.8e-6
print "fsSNE = ", fsSNE
To get comparable results with the USRP measurements, the number of samples was chosen so that the sensing time with SNE was approximately equal to the USRP sensing time.
In [4]:
ts = array([25.0e-3, 12.5e-3, 10.0e-3])
NsSNE = numpy.round(ts*fsSNE)
print "NsSNE =", NsSNE
The signal we attempted to detect was the IEEE Wireless Microphone simulation method, using the Soft Speaker profile.
IEEE Wireless Microphone simulation method simulates a transmission from a analog FM wireless microphone using an FM modulated sine wave:
$x(t) = A\cos\left[2\pi f_c t + 2 \pi \Delta f \int_0^t sin\left( 2\pi f_m u du\right) \right]$
$\Delta f = 15 \mathrm{kHz}$
$f_m = 3.9 \mathrm{kHz}$
The carrier frequency $f_c$ was 864 MHz for USRP (frequency for unlicensed operation of wireless microphones in Europe) and 850 MHz for SNE-ISMTV (frequency nearest to 864 MHz supported by the hardware).
Total attenuation in the cable, connectors and the 30 dB attenuator was measured using a R&S FSV spectrum analyzer. All equipment was left to warm-up for at least 2 hours before all measurements.
In [5]:
Image(filename='figures/power_measurement.png', width="100%")
Out[5]:
Figure above shows input power measured by the signal analyzer when the signal generator output is was set to -40 dBm.
In [6]:
Att = 71.53 - 40
print "Att =", Att, "dBm"
In [7]:
Pdmin = 0.9
Since in this controlled laboratory setup the noise power $P_{noise}$ remains constant for the chosen hardware configuration, the signal power directly relates to the signa-to-noise ratio (SNR)
$SNR = \frac{P_{in}}{P_{noise}}$
Here we assume that the coaxial cable and all equipment is well shielded, so that is no external interference. We also assume that the temperature and other factors affecting internal receiver noise do not change during the course of the experiment.
Before we can determine the probability of detection for a method we have to calcuate its detection threshold $\gamma_0$.
A probability of false alarm $P_{fa}$ was chosen.
In [8]:
Pfa = 0.1
Given $P_{fa}$ we can determine $\gamma_0$ from $\gamma$ complementary cumulative density function (CCDF) for the given detector with no input signal.
To measure the distribution of $\gamma$, the RF output of the generator was switched off and $N_p$ measurements of $\gamma$ taken for each of the chosen detectors.
In [9]:
Np = 1000
Given these measurements we approximate the CCDF and determine the value of $\gamma_0$ by interpolation.
In [10]:
def get_ccdf(x):
xs = array(x)
xs.sort()
N = float(len(xs))
P = arange(N)/N
return xs, P
def get_gamma0(gammaN):
gammaN, Pd = get_ccdf(gammaN)
gamma0 = interp(1. - Pfa, Pd, gammaN)
return gamma0
For example, energy detector using USRP with $f_s=1\mathrm{MHz}$ and $N_s=25000\mathrm{samples}$ yields the following value for the threshold:
In [11]:
gammaN = loadtxt("../simout-noisecomp-test/dat/sim_usrp_micsoft_fs1mhz_Ns25ks_ed_off.dat")
gamma0 = get_gamma0(gammaN)
gammaNs, Pd = get_ccdf(gammaN)
plot(gammaNs, 1.-Pd)
plot([gamma0, gamma0], [0, 1], 'b--')
plot([2e-5, 3e-5], [Pfa, Pfa], 'b--')
axis([2e-5, 3e-5, 0, 1])
title("CCDF")
xlabel("$\gamma_0$")
ylabel("$P_{fa}$")
grid()
In [12]:
print "gamma0 =", gamma0, "@ Pfa =", Pfa
In [13]:
def iterate_campaign(path):
for fn in glob.glob(path):
g = re.search("_m([0-9_]+)dbm\.dat$", fn)
if g:
Pg = -float(g.group(1).replace('_', '.'))
gamma = loadtxt(fn)
yield Pg, gamma
print "Pg =", array(sorted(Pg for Pg, gamma
in iterate_campaign("../simout-noisecomp-test/dat/sim_usrp_micsoft_fs1mhz_Ns25ks_ed_*.dat")))
If $P_g$ is the output power of the calibrated signal generator, then
$P_{in}$ = $P_g$ - $A$
Probability of detection $P_d$ is calculated for each $P_{in}$ and each method by counting the number of $\gamma$ measurements that were above the calculated threshold $\gamma_0$.
In [14]:
def get_campaign_g(path, gamma0):
Pg = []
Pd = []
for Pg0, gamma in iterate_campaign(path):
Pg.append(Pg0)
Pd.append(mean(gamma > gamma0))
Pg = array(Pg)
Pd = array(Pd)
Pga = Pg.argsort()
Pd = Pd[Pga]
Pg = Pg[Pga]
return Pg, Pd
def get_campaign(path, gamma0):
Pg, Pd = get_campaign_g(path, gamma0)
Pin = Pg - Att
return Pin, Pd
def get_Pinmin(path, Pdth=Pdmin):
gamma0 = None
for fn in glob.glob(path):
if fn.endswith("_off.dat"):
gammaN = loadtxt(fn)
gamma0 = get_gamma0(gammaN)
break
Pin, Pd = get_campaign(path, gamma0)
Pinmin = interp(Pdth, Pd, Pin)
return Pinmin
Given minimum $P_d$, $P_{dmin}$, the minimum signal power $P_{inmin}$ for the spectrum sensing method can be determined from measurements using interpolation.
For example, for energy detector using USRP:
In [15]:
Pinmin = get_Pinmin("../simout-noisecomp-test/dat/sim_usrp_micsoft_fs1mhz_Ns25ks_ed_*.dat")
print "Pinmin =", Pinmin, "@ Pdmin =", Pdmin
In [16]:
Pin, Pd = get_campaign("../simout-noisecomp-test/dat/sim_usrp_micsoft_fs1mhz_Ns25ks_ed_*.dat", gamma0)
plot(Pin, Pd, 'o-')
plot([Pinmin, Pinmin], [0, 1], 'b--')
x1, x2, y1, y2 = axis()
plot([x1, x2], [Pdmin, Pdmin], 'b--')
xlabel("$P_{in}$ [dBm]")
ylabel("$P_d$")
grid()
In [30]:
method_ps = set()
for path in glob.glob("../simout-noisecomp-test2/dat/sim_usrp_micsoft_*.dat"):
g = re.search("ks_(.*)_m", path)
if g:
method_ps.add(g.group(1))
method_ps = sorted(method_ps)
In [31]:
def get_batch_Pinmin_permethod(path, Pdth=Pdmin):
Pinmin_permethod = {}
for method_p in method_ps:
path2 = path.replace("*", "%s_*" % (method_p,))
Pinmin_permethod[method_p] = get_Pinmin(path2, Pdmin)
return Pinmin_permethod
In [44]:
from matplotlib import cm
def plot_comparisson(Pinmin_permethod, methods):
cmap = cm.get_cmap('gist_rainbow')
left = []
height = []
labels = []
colors = []
#methods = ['ed', 'scf', 'cav', 'cfn', 'mac', 'mme', 'eme', 'agm', 'met']
#methods = ['ed', 'cav', 'ccav', 'mac', 'cmac']
left1 = 0.
for n, method in enumerate(methods):
if method == 'ed':
method_ps = [ 'ismtv', 'ed']
elif method == 'scf':
method_ps = ["scf_Np%d" % scfNp for scfNp in [64, 128]]
else:
method_ps = [ "%s_l%d" % (method, l) for l in xrange(5, 25, 5) ]
for m, method_p in enumerate(method_ps):
try:
Pinmin = Pinmin_permethod[method_p]
except KeyError:
continue
left.append(left1)
height.append(Pinmin)
if method_p == 'ed':
label = 'ED (USRP)'
elif method_p == 'ismtv':
label = 'ED (ISMTV)'
else:
label = method_p.upper()
label = re.sub("_L(\d+)", r" $L=\1$", label)
label = re.sub("_NP(\d+)", r" $N'=\1$", label)
labels.append(label)
c = cmap(float(n)/(len(methods)-1) + m*.007)
# Color
colors.append(c)
# BW
#colors.append((.85,.85,.85))
left1 += 1.
left1 += .9
height = array(height)
left = array(left)
# print
#figure(figsize=(13, 5))
# screen
figure(figsize=(12, 5))
for x, y in zip(left, height):
text(x+.4, y-1.5, "%.1f" % y, horizontalalignment="center", rotation=80, backgroundcolor='w')
#o = 100
#bar(left, -height-o, color=colors, bottom=o)
bar(left, height, color=colors)
xticks(left+.4, labels, rotation=60, horizontalalignment='right')
ylabel("$P_{in}$ [dBm] @ $P_d = 0.9$")
#ylim(100, 118)
ylim(-130, -100)
xlim(-.5, left1-.5)
tight_layout()
grid()
# legend
yi = -129
for si, ei, label in [[0,2,"energy"],
[2,4,"cyclo-\nstationary"],
[4,20,"covariance"],
[20,32,"covariance eigenvalue"]]:
pass
#x1 = left[si]
#x2 = left[ei-1]+.8
#plot([x1,x2], [yi,yi], 'k')
#plot([x1,x1], [yi-.5,yi+.5], 'k')
#plot([x2,x2], [yi-.5,yi+.5], 'k')
#text((x1+x2)/2., yi+.8, label, horizontalalignment="center", backgroundcolor='w')
Lower detectable power is better
In [45]:
Pinmin_permethod_25ms = get_batch_Pinmin_permethod("../simout-noisecomp-test2/dat/sim_usrp_micsoft_fs1mhz_Ns25ks_*.dat")
Pinmin_permethod_25ms['ismtv'] = get_Pinmin("../measurements/pd/sneismtv/sim_sneismtv_micsoft_fs0mhz_Ns3ks_ed_n3676_*.dat")
methods = ['ed', 'scf', 'cav', 'cfn', 'mac', 'mme', 'eme', 'agm', 'met']
plot_comparisson(Pinmin_permethod_25ms, methods)
title("$t=25.0$ms (USRP @ $f_s=1$MHz, $N_s=25$ksamples, SNE-ISMTV @ $f_s=147$kHz, $N_s=3676$samples)");
#savefig("figures/pin_min_comparison_25ms.png", dpi=300)
#savefig("pin_min_comparison_25ms.eps")
None
In [46]:
Pinmin_permethod_25ms = get_batch_Pinmin_permethod("../simout-noisecomp-test2/dat/sim_usrp_micsoft_fs1mhz_Ns25ks_*.dat")
Pinmin_permethod_25ms['ismtv'] = get_Pinmin("../measurements/pd/sneismtv/sim_sneismtv_micsoft_fs0mhz_Ns3ks_ed_n3676_*.dat")
methods = ['ed', 'scf', 'ccav', 'ccfn', 'cmac', 'cmme', 'ceme', 'cagm', 'cmet']
plot_comparisson(Pinmin_permethod_25ms, methods)
title("$t=25.0$ms (USRP @ $f_s=1$MHz, $N_s=25$ksamples, SNE-ISMTV @ $f_s=147$kHz, $N_s=3676$samples, noise compensation)");
#savefig("figures/pin_min_comparison_25ms.png", dpi=300)
#savefig("pin_min_comparison_25ms.eps")
None
With a higher sampling rate, performance seems to decrease. Probably because less of the signal is covered with 25.000 samples.
In [47]:
Pinmin_permethod_13ms = get_batch_Pinmin_permethod("../simout-noisecomp-test2/dat/sim_usrp_micsoft_fs2mhz_Ns25ks_*.dat")
Pinmin_permethod_13ms['ismtv'] = get_Pinmin("../measurements/pd/sneismtv/sim_sneismtv_micsoft_fs0mhz_Ns3ks_ed_n1838_*.dat")
methods = ['ed', 'scf', 'cav', 'cfn', 'mac', 'mme', 'eme', 'agm', 'met']
plot_comparisson(Pinmin_permethod_13ms, methods)
title("$t=12.5$ms (USRP @ $f_s=2$MHz, $N_s=25$ksamples, SNE-ISMTV @ $f_s=147$kHz, $N_s=1838$samples)");
#savefig("figures/pin_min_comparison_13ms.png", dpi=300)
#savefig("pin_min_comparison_13ms.eps")
None
In [48]:
Pinmin_permethod_13ms = get_batch_Pinmin_permethod("../simout-noisecomp-test2/dat/sim_usrp_micsoft_fs2mhz_Ns25ks_*.dat")
Pinmin_permethod_13ms['ismtv'] = get_Pinmin("../measurements/pd/sneismtv/sim_sneismtv_micsoft_fs0mhz_Ns3ks_ed_n1838_*.dat")
methods = ['ed', 'scf', 'ccav', 'ccfn', 'cmac', 'cmme', 'ceme', 'cagm', 'cmet']
plot_comparisson(Pinmin_permethod_13ms, methods)
title("$t=12.5$ms (USRP @ $f_s=2$MHz, $N_s=25$ksamples, SNE-ISMTV @ $f_s=147$kHz, $N_s=1838$samples, noise compensation)");
#savefig("figures/pin_min_comparison_13ms.png", dpi=300)
#savefig("pin_min_comparison_13ms.eps")
None
In [51]:
Pinmin_permethod_10ms = get_batch_Pinmin_permethod("../simout-noisecomp-test2/dat/sim_usrp_micsoft_fs10mhz_Ns100ks_*.dat")
Pinmin_permethod_10ms['ismtv'] = get_Pinmin("../measurements/pd/sneismtv/sim_sneismtv_micsoft_fs0mhz_Ns3ks_ed_n1471_*.dat")
methods = ['ed', 'scf', 'cav', 'cfn', 'mac', 'mme', 'eme', 'agm', 'met']
plot_comparisson(Pinmin_permethod_10ms, methods)
title("$t=10.0$ms (USRP @ $f_s=10$MHz, $N_s=100$ksamples, SNE-ISMTV @ $f_s=147$kHz, $N_s=1471$samples)");
#savefig("figures/pin_min_comparison_10ms.png", dpi=300)
#savefig("pin_min_comparison_10ms.eps")
None
In [52]:
Pinmin_permethod_10ms = get_batch_Pinmin_permethod("../simout-noisecomp-test2/dat/sim_usrp_micsoft_fs10mhz_Ns100ks_*.dat")
Pinmin_permethod_10ms['ismtv'] = get_Pinmin("../measurements/pd/sneismtv/sim_sneismtv_micsoft_fs0mhz_Ns3ks_ed_n1471_*.dat")
methods = ['ed', 'scf', 'ccav', 'ccfn', 'cmac', 'cmme', 'ceme', 'cagm', 'cmet']
plot_comparisson(Pinmin_permethod_10ms, methods)
title("$t=10.0$ms (USRP @ $f_s=10$MHz, $N_s=100$ksamples, SNE-ISMTV @ $f_s=147$kHz, $N_s=1471$samples, noise compensation)");
#savefig("figures/pin_min_comparison_10ms.png", dpi=300)
#savefig("pin_min_comparison_10ms.eps")
None
In [ ]:
In [ ]: